Material didáctico de: Modelos Supervisados y Evaluación de Modelos.
Complementarios:
En este ejercicio vamos a seguir los pasos del ciclo de vida de un proyecto de minería de datos, para el caso de un algoritmo de clasificación y más concretamente un árbol de decisión. Primero y a modo de ejemplo sencillo lo haremos con el archivo titanic.csv, que se encuentra adjunto en el aula. Este archivo contiene un registro por cada pasajero que viajaba en el Titanic. En las variables se caracteriza si era hombre o mujer, adulto o menor (niño), en qué categoría viajaba o si era miembro de la tripulación. Se mostrará un ejemplo sencillo de solución con estos datos pero los alumnos deberéis responder a las preguntas de la rúbrica para otro conjunto: German Credit. Para este conjunto, tomaréis como referencia la variable “default” que indica el impago de créditos.
Objetivos:
A continuación, se plantean los puntos a realizar en la PEC 3 y, tomando como ejemplo el conjunto de datos de Titanic, se obtendrán, a modo de ejemplo, algunos resultados que pretender servir a modo de inspiración para los estudiantes. Los estudiantes deberán utilizar el conjunto de datos de “German Credit Data” que se pueden conseguir en este enlace: https://www.kaggle.com/shravan3273/credit-approval
Revisión de los datos, extracción visual de información y preparación de los datos
Carga de los datos:
data<-read.csv("./titanic.csv",header=T,sep=",")
attach(data)
Empezaremos haciendo un breve análisis de los datos ya que nos interesa tener una idea general de los datos que disponemos.
Primero calcularemos las dimensiones de nuestra base de datos y analizaremos qué tipos de atributos tenemos.
Para empezar, calculamos las dimensiones de la base de datos mediante la función dim(). Obtenemos que disponemos de 2201 registros o pasajeros (filas) y 4 variables (columnas).
dim(data)
## [1] 2201 4
¿Cuáles son esas variables? Gracias a la función str() sabemos que las cuatro variables son categóricas o discretas, es decir, toman valores en un conjunto finito. La variable CLASS hace referencia a la clase en la que viajaban los pasajeros (1ª, 2ª, 3ª o crew), AGE determina si era adulto o niño (Adulto o Menor), la variable SEX si era hombre o mujer (Hombre o Mujer) y la última variable (SURVIVED) informa si el pasajero murió o sobrevivió en el accidente (Muere o Sobrevive).
str(data)
## 'data.frame': 2201 obs. of 4 variables:
## $ CLASS : chr "1a" "1a" "1a" "1a" ...
## $ AGE : chr "Adulto" "Adulto" "Adulto" "Adulto" ...
## $ SEX : chr "Hombre" "Hombre" "Hombre" "Hombre" ...
## $ SURVIVED: chr "Sobrevive" "Sobrevive" "Sobrevive" "Sobrevive" ...
Es de gran interés saber si tenemos muchos valores nulos (campos vacíos) y la distribución de valores por variables. Es por ello recomendable empezar el análisis con una visión general de las variables. Mostraremos para cada atributo la cantidad de valores perdidos mediante la función summary.
summary(data)
## CLASS AGE SEX SURVIVED
## Length:2201 Length:2201 Length:2201 Length:2201
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
Como parte de la preparación de los datos, miraremos si hay valores missing.
missing <- data[is.na(data),]
dim(missing)
## [1] 0 4
Observamos fácilmente que no hay valores missing y, por tanto, no deberemos preparar los datos en este sentido. En caso de haberlos, habría que tomar decisiones para tratar los datos adecuadamente.
Disponemos por tanto de un data frame formado por cuatro variables categóricas sin valores nulos.
Para un conocimiento mayor sobre los datos, tenemos a nuestro alcance unas herramientas muy valiosas: las herramientas de visualización. Para dichas visualizaciones, haremos uso de los paquetes ggplot2, gridExtra y grid de R.
if(!require(ggplot2)){
install.packages('ggplot2', repos='http://cran.us.r-project.org')
library(ggplot2)
}
## Loading required package: ggplot2
if(!require(ggpubr)){
install.packages('ggpubr', repos='http://cran.us.r-project.org')
library(ggpubr)
}
## Loading required package: ggpubr
if(!require(grid)){
install.packages('grid', repos='http://cran.us.r-project.org')
library(grid)
}
## Loading required package: grid
if(!require(gridExtra)){
install.packages('gridExtra', repos='http://cran.us.r-project.org')
library(gridExtra)
}
## Loading required package: gridExtra
if(!require(C50)){
install.packages('C50', repos='http://cran.us.r-project.org')
library(C50)
}
## Loading required package: C50
Siempre es importante analizar los datos que tenemos ya que las conclusiones dependerán de las características de la muestra.
grid.newpage()
plotbyClass<-ggplot(data,aes(CLASS))+geom_bar() +labs(x="Class", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("blue","#008000"))+ggtitle("Class")
plotbyAge<-ggplot(data,aes(AGE))+geom_bar() +labs(x="Age", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("blue","#008000"))+ggtitle("Age")
plotbySex<-ggplot(data,aes(SEX))+geom_bar() +labs(x="Sex", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("blue","#008000"))+ggtitle("Sex")
plotbySurvived<-ggplot(data,aes(SURVIVED))+geom_bar() +labs(x="Survived", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("blue","#008000"))+ggtitle("SURVIVED")
grid.arrange(plotbyClass,plotbyAge,plotbySex,plotbySurvived,ncol=2)
Claramente vemos cómo es la muestra analizando la distribución de las variables disponibles. De cara a los informes, es mucho más interesante esta información que la obtenida en summary, que se puede usar para complementar.
Nos interesa describir la relación entre la supervivencia y cada uno de las variables mencionadas anteriormente. Para ello, por un lado graficaremos mediante diagramas de barras la cantidad de muertos y supervivientes según la clase en la que viajaban, la edad o el sexo. Por otro lado, para obtener los datos que estamos graficando utilizaremos el comando table para dos variables que nos proporciona una tabla de contingencia.
grid.newpage()
plotbyClass<-ggplot(data,aes(CLASS,fill=SURVIVED))+geom_bar() +labs(x="Class", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("black","#008000"))+ggtitle("Survived by Class")
plotbyAge<-ggplot(data,aes(AGE,fill=SURVIVED))+geom_bar() +labs(x="Age", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("black","#008000"))+ggtitle("Survived by Age")
plotbySex<-ggplot(data,aes(SEX,fill=SURVIVED))+geom_bar() +labs(x="Sex", y="Passengers")+ guides(fill=guide_legend(title=""))+ scale_fill_manual(values=c("black","#008000"))+ggtitle("Survived by Sex")
grid.arrange(plotbyClass,plotbyAge,plotbySex,ncol=2)
De estos gráficos obtenemos información muy valiosa que complementamos con las tablas de contingencia (listadas abajo). Por un lado, la cantidad de pasajeros que sobrevivieron es similar en hombres y mujeres (hombres: 367 y mujeres 344). No, en cambio, si tenemos en cuenta el porcentaje respecto a su sexo. Es decir, pese a que la cantidad de mujeres y hombres que sobrevivieron es pareja, viajaban más hombres que mujeres (470 mujeres y 1731 hombres), por lo tanto, la tasa de muerte en hombres es muchísimo mayor (el 78,79% de los hombres murieron mientras que en mujeres ese porcentaje baja a 26,8%).
En cuanto a la clase en la que viajaban, los pasajeros que viajaban en primera clase fueron los únicos que el porcentaje de supervivencia era mayor que el de mortalidad. El 62,46% de los viajeros de primera clase sobrevivió, el 41,4% de los que viajaban en segunda clase mientras que de los viajeros de tercera y de la tripulación solo sobrevivieron un 25,21% y 23,95% respectivamente. Para finalizar, destacamos que la presencia de pasajeros adultos era mucho mayor que la de los niños (2092 frente a 109) y que la tasa de supervivencia en niños fue mucho mayor (52,29% frente a 31,26%), no podemos obviar, en cambio, que los únicos niños que murieron fueron todos pasajeros de tercera clase (52 niños).
tabla_SST <- table(SEX, SURVIVED)
tabla_SST
## SURVIVED
## SEX Muere Sobrevive
## Hombre 1364 367
## Mujer 126 344
prop.table(tabla_SST, margin = 1)
## SURVIVED
## SEX Muere Sobrevive
## Hombre 0.7879838 0.2120162
## Mujer 0.2680851 0.7319149
tabla_SCT <- table(CLASS,SURVIVED)
tabla_SCT
## SURVIVED
## CLASS Muere Sobrevive
## 1a 122 203
## 2a 167 118
## 3a 528 178
## crew 673 212
prop.table(tabla_SCT, margin = 1)
## SURVIVED
## CLASS Muere Sobrevive
## 1a 0.3753846 0.6246154
## 2a 0.5859649 0.4140351
## 3a 0.7478754 0.2521246
## crew 0.7604520 0.2395480
tabla_SAT <- table(AGE,SURVIVED)
tabla_SAT
## SURVIVED
## AGE Muere Sobrevive
## Adulto 1438 654
## Menor 52 57
prop.table(tabla_SAT, margin = 1)
## SURVIVED
## AGE Muere Sobrevive
## Adulto 0.6873805 0.3126195
## Menor 0.4770642 0.5229358
tabla_SAT.byClass <- table(AGE,SURVIVED,CLASS)
tabla_SAT.byClass
## , , CLASS = 1a
##
## SURVIVED
## AGE Muere Sobrevive
## Adulto 122 197
## Menor 0 6
##
## , , CLASS = 2a
##
## SURVIVED
## AGE Muere Sobrevive
## Adulto 167 94
## Menor 0 24
##
## , , CLASS = 3a
##
## SURVIVED
## AGE Muere Sobrevive
## Adulto 476 151
## Menor 52 27
##
## , , CLASS = crew
##
## SURVIVED
## AGE Muere Sobrevive
## Adulto 673 212
## Menor 0 0
Los resultados anteriores muestran los datos de forma descriptiva, podemos añadir algún test estadístico para validar el grado de significancia de la relación. La librería “DescTools” nos permite instalarlo fácilmente.
if(!require(DescTools)){
install.packages('DescTools', repos='http://cran.us.r-project.org')
library(DescTools)
}
## Loading required package: DescTools
Phi(tabla_SST)
## [1] 0.4556048
CramerV(tabla_SST)
## [1] 0.4556048
Phi(tabla_SAT)
## [1] 0.09757511
CramerV(tabla_SAT)
## [1] 0.09757511
Phi(tabla_SCT)
## [1] 0.2941201
CramerV(tabla_SCT)
## [1] 0.2941201
Valores de la V de Cramér (https://en.wikipedia.org/wiki/Cramér%27s_V) y Phi (https://en.wikipedia.org/wiki/Phi_coefficient) entre 0.1 y 0.3 nos indican que la asociación estadística es baja, y entre 0.3 y 0.5 se puede considerar una asociación media. Finalmente, si los valores fueran superiores a 0.5 (no es el caso), la asociación estadística entre las variables sería alta. Como se puede apreciar, los valores de Phi y V coinciden. Esto ocurre en el contexto de analizar tablas de contingencia 2x2.
Una alternativa interesante a las barras de diagramas, es el plot de las tablas de contingencia. Obtenemos la misma información pero para algunos receptores puede resultar más visual.
par(mfrow=c(2,2))
plot(tabla_SCT, col = c("black","#008000"), main = "SURVIVED vs. CLASS")
plot(tabla_SAT, col = c("black","#008000"), main = "SURVIVED vs. AGE")
plot(tabla_SST, col = c("black","#008000"), main = "SURVIVED vs. SEX")
Nuestro objetivo es crear un árbol de decisión que permita analizar qué tipo de pasajero del Titanic tenía probabilidades de sobrevivir o no. Por lo tanto, la variable por la que clasificaremos es el campo de si el pasajero sobrevivió o no. De todas maneras, al imprimir las primeras (con head) y últimas 10 (con tail) filas nos damos cuenta de que los datos están ordenados.
head(data,10)
## CLASS AGE SEX SURVIVED
## 1 1a Adulto Hombre Sobrevive
## 2 1a Adulto Hombre Sobrevive
## 3 1a Adulto Hombre Sobrevive
## 4 1a Adulto Hombre Sobrevive
## 5 1a Adulto Hombre Sobrevive
## 6 1a Adulto Hombre Sobrevive
## 7 1a Adulto Hombre Sobrevive
## 8 1a Adulto Hombre Sobrevive
## 9 1a Adulto Hombre Sobrevive
## 10 1a Adulto Hombre Sobrevive
tail(data,10)
## CLASS AGE SEX SURVIVED
## 2192 crew Adulto Mujer Sobrevive
## 2193 crew Adulto Mujer Sobrevive
## 2194 crew Adulto Mujer Sobrevive
## 2195 crew Adulto Mujer Sobrevive
## 2196 crew Adulto Mujer Sobrevive
## 2197 crew Adulto Mujer Sobrevive
## 2198 crew Adulto Mujer Sobrevive
## 2199 crew Adulto Mujer Muere
## 2200 crew Adulto Mujer Muere
## 2201 crew Adulto Mujer Muere
Nos interesa “desordenarlos”. Guardaremos los datos con el nuevo nombre como “data_random”.
set.seed(1)
data_random <- data[sample(nrow(data)),]
Para la futura evaluación del árbol de decisión, es necesario dividir el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba. El conjunto de entrenamiento es el subconjunto del conjunto original de datos utilizado para construir un primer modelo; y el conjunto de prueba, el subconjunto del conjunto original de datos utilizado para evaluar la calidad del modelo.
Lo más correcto será utilizar un conjunto de datos diferente del que utilizamos para construir el árbol, es decir, un conjunto diferente del de entrenamiento. No hay ninguna proporción fijada con respecto al número relativo de componentes de cada subconjunto, pero la más utilizada acostumbra a ser 2/3 para el conjunto de entrenamiento y 1/3, para el conjunto de prueba.
La variable por la que clasificaremos es el campo de si el pasajero sobrevivió o no, que está en la cuarta columna. De esta forma, tendremos un conjunto de datos para el entrenamiento y uno para la validación
set.seed(666)
y <- data_random[,4]
X <- data_random[,1:3]
De forma dinámica podemos definir una forma de separar los datos en función de un parámetro, en este caso del “split_prop”. Definimos un parámetro que controla el split de forma dinámica en el test.
split_prop <- 3
max_split<-floor(nrow(X)/split_prop)
tr_limit <- nrow(X)-max_split
ts_limit <- nrow(X)-max_split+1
trainX <- X[1:tr_limit,]
trainy <- y[1:tr_limit]
testX <- X[ts_limit+1:nrow(X),]
testy <- y[ts_limit+1:nrow(X)]
En la segunda opción podemos crear directamente un rango utilizando el mismo parámetro anterior.
split_prop <- 3
indexes = sample(1:nrow(data), size=floor(((split_prop-1)/split_prop)*nrow(data)))
trainX<-X[indexes,]
trainy<-y[indexes]
testX<-X[-indexes,]
testy<-y[-indexes]
Después de una extracción aleatoria de casos es altamente recomendable efectuar un análisis de datos mínimo para asegurarnos de no obtener clasificadores sesgados por los valores que contiene cada muestra. En este caso, verificaremos que la proporción del supervivientes es más o menos constante en los dos conjuntos:
summary(trainX);
## CLASS AGE SEX
## Length:1467 Length:1467 Length:1467
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
summary(trainy)
## Length Class Mode
## 1467 character character
summary(testX)
## CLASS AGE SEX
## Length:734 Length:734 Length:734
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
summary(testy)
## Length Class Mode
## 734 character character
Verificamos fácilmente que no hay diferencias graves que puedan sesgar las conclusiones.
Se crea el árbol de decisión usando los datos de entrenamiento (no hay que olvidar que la variable outcome es de tipo factor):
trainy = as.factor(trainy)
model <- C50::C5.0(trainX, trainy,rules=TRUE )
summary(model)
##
## Call:
## C5.0.default(x = trainX, y = trainy, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Sat Dec 11 10:50:00 2021
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 1467 cases (4 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (1169/255, lift 1.2)
## SEX = Hombre
## -> class Muere [0.781]
##
## Rule 2: (20, lift 2.9)
## CLASS in {2a, 1a}
## AGE = Menor
## -> class Sobrevive [0.955]
##
## Rule 3: (298/75, lift 2.3)
## SEX = Mujer
## -> class Sobrevive [0.747]
##
## Default class: Muere
##
##
## Evaluation on training data (1467 cases):
##
## Rules
## ----------------
## No Errors
##
## 3 317(21.6%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 914 75 (a): class Muere
## 242 236 (b): class Sobrevive
##
##
## Attribute usage:
##
## 100.00% SEX
## 1.36% CLASS
## 1.36% AGE
##
##
## Time: 0.0 secs
Errors muestra el número y porcentaje de casos mal clasificados en el subconjunto de entrenamiento. El árbol obtenido clasifica erróneamente 317 de los 1467 casos dados, una tasa de error del 21.6%.
A partir del árbol de decisión de dos hojas que hemos modelado, se pueden extraer las siguientes reglas de decisión (gracias a rules=TRUE podemos imprimir las reglas directamente):
SEX = “Hombre” → Muere. Validez: 78,1%
CLASS “1ª”, “2ª” y AGE = “Menor” → Sobrevive. Validez: 95,5%
SEX = “Mujer” → Sobrevive. Validez: 74,7%
Por tanto, podemos concluir que el conocimiento extraído y cruzado con el análisis visual se resume en “las mujeres y los niños primero a excepción de que fueras de 3ª clase”.
A continuación, mostramos el árbol obtenido.
model <- C50::C5.0(trainX, trainy)
plot(model)
Una vez tenemos el modelo, podemos comprobar su calidad prediciendo la clase para los datos de prueba que nos hemos reservado al principio.
predicted_model <- predict( model, testX, type="class" )
print(sprintf("La precisión del árbol es: %.4f %%",100*sum(predicted_model == testy) / length(predicted_model)))
## [1] "La precisión del árbol es: 78.2016 %"
Cuando hay pocas clases, la calidad de la predicción se puede analizar mediante una matriz de confusión que identifica los tipos de errores cometidos.
mat_conf<-table(testy,Predicted=predicted_model)
mat_conf
## Predicted
## testy Muere Sobrevive
## Muere 450 51
## Sobrevive 109 124
Otra manera de calcular el porcentaje de registros correctamente clasificados usando la matriz de confusión:
porcentaje_correct<-100 * sum(diag(mat_conf)) / sum(mat_conf)
print(sprintf("El %% de registros correctamente clasificados es: %.4f %%",porcentaje_correct))
## [1] "El % de registros correctamente clasificados es: 78.2016 %"
Además, tenemos a nuestra disposición el paquete gmodels para obtener información más completa:
if(!require(gmodels)){
install.packages('gmodels', repos='http://cran.us.r-project.org')
library(gmodels)
}
## Loading required package: gmodels
## Registered S3 method overwritten by 'gdata':
## method from
## reorder.factor DescTools
CrossTable(testy, predicted_model,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 734
##
##
## | Prediction
## Reality | Muere | Sobrevive | Row Total |
## -------------|-----------|-----------|-----------|
## Muere | 450 | 51 | 501 |
## | 0.613 | 0.069 | |
## -------------|-----------|-----------|-----------|
## Sobrevive | 109 | 124 | 233 |
## | 0.149 | 0.169 | |
## -------------|-----------|-----------|-----------|
## Column Total | 559 | 175 | 734 |
## -------------|-----------|-----------|-----------|
##
##
En este apartado buscaremos probar con las variaciones que nos ofrece el paquete C5.0 para analizar cómo afectan a la creación de los árboles generados. Existen muchas posibles variaciones con otras funciones que podéis investigar. La idea es seguir con el enfoque de árboles de decisión explorando posibles opciones. Una vez tengamos un método alternativo, debemos analizar cómo se modifica el árbol y cómo afecta a la capacidad predictiva en el conjunto de test.
A continuación, utilizamos otro enfoque para comparar los resultados: incorpora como novedad “adaptative boosting”, basado en el trabajo Rob Schapire and Yoav Freund (1999). La idea de esta técnica es generar varios clasificadores, con sus correspondientes arboles de decisión y su ser de reglas. Cuando un nuevo caso va a ser clasificado, cada clasificador vota cual es la clase predicha. Los votos son sumados y determina la clase final.
modelo2 <- C50::C5.0(trainX, trainy, trials = 10)
plot(modelo2)
En este caso, dada la simplicidad del conjunto de ejemplo, no se aprecian diferencias, pero aparecerán en datos de mayor complejidad y modificando el parámetro “trials” se puede intentar mejorar los resultados.
Vemos a continuación cómo son las predicciones del nuevo árbol:
predicted_model2 <- predict( modelo2, testX, type="class" )
print(sprintf("La precisión del árbol es: %.4f %%",100*sum(predicted_model2 == testy) / length(predicted_model2)))
## [1] "La precisión del árbol es: 79.5640 %"
Observamos como se modifica levemente la precisión del modelo a mejor.
mat_conf<-table(testy,Predicted=predicted_model2)
mat_conf
## Predicted
## testy Muere Sobrevive
## Muere 491 10
## Sobrevive 140 93
Otra manera de calcular el porcentaje de registros correctamente clasificados usando la matriz de confusión:
porcentaje_correct<-100 * sum(diag(mat_conf)) / sum(mat_conf)
print(sprintf("El %% de registros correctamente clasificados es: %.4f %%",porcentaje_correct))
## [1] "El % de registros correctamente clasificados es: 79.5640 %"
El algoritmo C5.0 incorpora algunas opciones para ver la importancia de las variables (ver documentación para los detalles entre los dos métodos):
importancia_usage <- C50::C5imp(modelo2, metric = "usage")
importancia_splits <- C50::C5imp(modelo2, metric = "splits")
importancia_usage
## Overall
## CLASS 100.00
## SEX 100.00
## AGE 91.55
importancia_splits
## Overall
## CLASS 36.84211
## SEX 36.84211
## AGE 26.31579
Curiosamente y aunque el conjunto de datos es muy sencillo, se aprecian diferencias en los métodos de importancia de las variables. Se recomienda en vuestro ejercicio mejorar la visualización de los resultados con la función ggplo2 o similar.
Para el conjunto de datos German Credit, los alumnos deben completar aquí la solución a la PEC3 que consiste de los siguientes apartados. Notad que se detalla el contenido necesario para cada apartado en la Sección 4 (Rúbrica).
Carga de los datos:
datos<-read.csv(file = "./credit.csv",header=T,sep=",")
Cargar librerias:
install.packages("fBasics")
## package 'fBasics' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\eduar\AppData\Local\Temp\RtmpmI0AcG\downloaded_packages
library(fBasics) #summary statistics
library(plotly) #data visulization
library(dplyr) #count function
library(C50) #C50 Decision Tree algorithm
library(gmodels) #CrossTable()
library(caret) #Confusion Matrix
library(rmarkdown)
Lo primero que debemos hacer es analizar las variables que componen el fichero, para ello obtendremos su dimensión, el número de filas, su estructura y sus estadísticas básicas.
#Dimensión
dim(datos)
## [1] 1000 21
#Filas
filas=dim(datos)[1]
Podemos ver que tenemos 1000 entradas distintas y tenemos 21 características que definen el conjunto de datos
#Estructura
str(datos)
## 'data.frame': 1000 obs. of 21 variables:
## $ checking_balance : chr "< 0 DM" "1 - 200 DM" "unknown" "< 0 DM" ...
## $ months_loan_duration: int 6 48 12 42 24 36 24 36 12 30 ...
## $ credit_history : chr "critical" "repaid" "critical" "repaid" ...
## $ purpose : chr "radio/tv" "radio/tv" "education" "furniture" ...
## $ amount : int 1169 5951 2096 7882 4870 9055 2835 6948 3059 5234 ...
## $ savings_balance : chr "unknown" "< 100 DM" "< 100 DM" "< 100 DM" ...
## $ employment_length : chr "> 7 yrs" "1 - 4 yrs" "4 - 7 yrs" "4 - 7 yrs" ...
## $ installment_rate : int 4 2 2 2 3 2 3 2 2 4 ...
## $ personal_status : chr "single male" "female" "single male" "single male" ...
## $ other_debtors : chr "none" "none" "none" "guarantor" ...
## $ residence_history : int 4 2 3 4 4 4 4 2 4 2 ...
## $ property : chr "real estate" "real estate" "real estate" "building society savings" ...
## $ age : int 67 22 49 45 53 35 53 35 61 28 ...
## $ installment_plan : chr "none" "none" "none" "none" ...
## $ housing : chr "own" "own" "own" "for free" ...
## $ existing_credits : int 2 1 1 1 2 1 1 1 1 2 ...
## $ default : int 1 2 1 1 2 1 1 1 1 2 ...
## $ dependents : int 1 1 2 2 2 2 1 1 1 1 ...
## $ telephone : chr "yes" "none" "none" "none" ...
## $ foreign_worker : chr "yes" "yes" "yes" "yes" ...
## $ job : chr "skilled employee" "skilled employee" "unskilled resident" "skilled employee" ...
Viendo la estructura vemos que tenemos 7 variables de tipo numérico y el resto son de tipo categórico.
#Estadísticas básicas
summary(datos)
## checking_balance months_loan_duration credit_history purpose amount savings_balance employment_length installment_rate personal_status
## Length:1000 Min. : 4.0 Length:1000 Length:1000 Min. : 250 Length:1000 Length:1000 Min. :1.000 Length:1000
## Class :character 1st Qu.:12.0 Class :character Class :character 1st Qu.: 1366 Class :character Class :character 1st Qu.:2.000 Class :character
## Mode :character Median :18.0 Mode :character Mode :character Median : 2320 Mode :character Mode :character Median :3.000 Mode :character
## Mean :20.9 Mean : 3271 Mean :2.973
## 3rd Qu.:24.0 3rd Qu.: 3972 3rd Qu.:4.000
## Max. :72.0 Max. :18424 Max. :4.000
## other_debtors residence_history property age installment_plan housing existing_credits default dependents
## Length:1000 Min. :1.000 Length:1000 Min. :19.00 Length:1000 Length:1000 Min. :1.000 Min. :1.0 Min. :1.000
## Class :character 1st Qu.:2.000 Class :character 1st Qu.:27.00 Class :character Class :character 1st Qu.:1.000 1st Qu.:1.0 1st Qu.:1.000
## Mode :character Median :3.000 Mode :character Median :33.00 Mode :character Mode :character Median :1.000 Median :1.0 Median :1.000
## Mean :2.845 Mean :35.55 Mean :1.407 Mean :1.3 Mean :1.155
## 3rd Qu.:4.000 3rd Qu.:42.00 3rd Qu.:2.000 3rd Qu.:2.0 3rd Qu.:1.000
## Max. :4.000 Max. :75.00 Max. :4.000 Max. :2.0 Max. :2.000
## telephone foreign_worker job
## Length:1000 Length:1000 Length:1000
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
Y gracias a las estadísticas básicas, podemos ver (sobre todo en las variables de tipo numérica) entre que rangos nos encontramos:
–> checking_balance: Estado de la cuenta corriente existente. Los valores que esta variable puede tomar son: <0 DM, 1 - 200 DM, unknown o > 200 DM.
–> months_loan_duration: Duración del crédito en meses.
–> credit_history: Historial de crédito. Los valores que esta variable puede tomar son: repaid, fully repaid this bank, fully repaid, delayed y critical.
–> purpose: Propósito. Los valores que esta variable puede tomar son: automóvil (nuevo), automóvil (usado), mobiliario / equipo, radio / televisión, electrodomésticos, reparaciones, educación, reciclaje, empresa u otros.
–> amount: Importe del crédito.
–> savings_balance: Cantidad en la cuenta de ahorro. Los valores que esta variable puede tomar son: < 100 DM, 101 - 500 DM, 501 - 1000 DM o unknown.
–> employment_length: Duración del empleo actual: > 7 yrs, 0 - 1 yrs, 1 - 4 yrs, 4 - 7 yrs y unemployed.
–> installment_rate: Tasa de cuotas en porcentaje de la renta disponible.
–> personal_status: Estado personal y sexo. Los valores que esta variable puede tomar son: hombre: divorciado / separado, mujer: divorciado / separado / casado, hombre: soltero, hombre: casado / viudo y mujer: soltero.
–> other_debtors: Otros deudores. Los valores que esta variable puede tomar son: none, co-applicant y guarantor.
–> residence_history: Residencia actual desde.
–> property: Propiedades o seguros. Los valores que esta variable puede tomar son: building society savings, other, real estate y unknown/none.
–> age: Edad en años.
–> installment_plan: Otros planes de pago. Los valores que esta variable puede tomar son: bank, none y stores.
–> housing: Tipo de vivienda. Los valores que esta variable puede tomar son: for free, own y rent.
–> existing_credits: Número de créditos existentes en este banco.
–> default: ¿Ha pagado el crédito?. Los valores que esta variable puede tomar son: Yes o No.
–> dependents: Número de personas que están obligadas a proporcionar mantenimiento para el crédito.
–> telephone: Teléfono. Los valores que esta variable puede tomar son: Yes o No.
–> foreign_worker: trabajador en el extranjero. Los valores que esta variable puede tomar son: Yes o No.
–> job: Tipo de trabajo mangement. Los valores que esta variable puede tomar son: self-employed, skilled employee, unemployed non-resident y unskilled resident.
Lo primero que vamos a hacer es normalizar el campo default:
datos$default[datos$default == 1] <- "NO MOROSO"
datos$default[datos$default == 2] <- "MOROSO"
datos$default <- as.factor(datos$default)
A continuación, se crearán algunas tablas y gráficas útiles que ofrecen estadísticas resumidas de las características del préstamo, como cuenta corriente, ahorros, duración y monto en función del incumplimiento (checking, savings, duration, amount, and default).
#Gráfica
datos %>% count(default, checking_balance) %>% plot_ly(x = ~default, y = ~n, color = ~checking_balance, type = "bar") %>% layout(title = "CANTIDAD DE SALDO EN LA CUENTA VS CREDITO SALDADO", xaxis = list(title = "CANTIDAD DE SALDO EN LA CUENTA"), yaxis = list(title = "CANTIDAD"))
#Estadisticas
knitr::kable(data.frame(table(datos$checking_balance)), Caption = "Checking Balance", col.names = c("Cuenta corriente", "Frecuencia"))
| Cuenta corriente | Frecuencia |
|---|---|
| < 0 DM | 274 |
| > 200 DM | 63 |
| 1 - 200 DM | 269 |
| unknown | 394 |
Como se puede observar, las cuentas desconocidas son la inmensa mayoría, pero las cuentas de entre 1 - 200 DM son las que suelen pagar más los créditos.
#Gráfica
datos %>% count(default, savings_balance) %>% plot_ly(x = ~default, y = ~n, color = ~savings_balance, type = "bar") %>% layout(title = "CANTIDAD DE AHORRO VS CREDITO SALDADO", xaxis = list(title = "CANTIDAD DE AHORRO"), yaxis = list(title = "CANTIDAD"))
#Estadisticas
knitr::kable(data.frame(table(datos$savings_balance)), Caption = "savings_balance", col.names = c("Saldo ahorrado", "Frecuencia"))
| Saldo ahorrado | Frecuencia |
|---|---|
| < 100 DM | 603 |
| > 1000 DM | 48 |
| 101 - 500 DM | 103 |
| 501 - 1000 DM | 63 |
| unknown | 183 |
Como se puede observar, las personas suelen tener un saldo ahorrado < 100 DM, y tienen un índice menor de morosidad.
#Gráfica
plot_ly(datos, y = ~months_loan_duration, color = ~default, type = "box") %>% layout(title = "DURACIÓN CRÉDITO VS CREDITO SALDADO")
#Estadisticas
knitr::kable(basicStats(datos$months_loan_duration), Caption = "Duración", digits = 2, col.names = c("Valores"))
| Valores | |
|---|---|
| nobs | 1000.00 |
| NAs | 0.00 |
| Minimum | 4.00 |
| Maximum | 72.00 |
| 1. Quartile | 12.00 |
| 3. Quartile | 24.00 |
| Mean | 20.90 |
| Median | 18.00 |
| Sum | 20903.00 |
| SE Mean | 0.38 |
| LCL Mean | 20.15 |
| UCL Mean | 21.65 |
| Variance | 145.42 |
| Stdev | 12.06 |
| Skewness | 1.09 |
| Kurtosis | 0.90 |
Como se puede observar, los créditos más cortos se suelen pagar más que los largos Además, la tabla de estadísticas da una información bastante completa sobre los cuartiles, la media, varianza…
#Gráfica
plot_ly(datos, y = ~amount, color = ~default, type = "box") %>% layout(title = "CANTIDAD CRÉDITO VS CREDITO SALDADO")
#Estadisticas
knitr::kable(basicStats(datos$amount), Caption = "Importe del crédito", digits = 2, col.names = c("Valores"))
| Valores | |
|---|---|
| nobs | 1000.00 |
| NAs | 0.00 |
| Minimum | 250.00 |
| Maximum | 18424.00 |
| 1. Quartile | 1365.50 |
| 3. Quartile | 3972.25 |
| Mean | 3271.26 |
| Median | 2319.50 |
| Sum | 3271258.00 |
| SE Mean | 89.26 |
| LCL Mean | 3096.09 |
| UCL Mean | 3446.42 |
| Variance | 7967843.47 |
| Stdev | 2822.74 |
| Skewness | 1.94 |
| Kurtosis | 4.25 |
Como se puede observar, los créditos más pequeños se suelen pagar más que los más grandes Además, la tabla de estadísticas da una información bastante completa sobre los cuartiles, la media, varianza…
A primera vista podemos observar que los créditos mas grandes y la cantidad de tiempo de pago mas larga son los que menos se suelen pagar. Además, las personas con una cantidad de dinero ahorrado o en la cuenta menor son los más morosos.
Ahora seguiremos explorando el resto de los campos, visualizándolos y comparándolos con el campo default.
Los siguientes campos para comparar son los relativos a la personas y su trabajo:
#Gráfica EDAD
plot_ly(datos, y = ~age, color = ~default, type = "box") %>% layout(title = "EDAD VS CREDITO SALDADO")
#Gráfica tipo Trabajo
datos %>% count(default, job) %>% plot_ly(x = ~default, y = ~n, color = ~job, type = "bar") %>% layout(title = "TIPO DE TRABAJO VS CREDITO SALDADO", xaxis = list(title = "TIPO DE TRABAJO"), yaxis = list(title = "CANTIDAD"))
#Gráfica duración del Trabajo
datos %>% count(default, employment_length) %>% plot_ly(x = ~default, y = ~n, color = ~employment_length, type = "bar") %>% layout(title = "DURACIÓN TRABAJO VS CREDITO SALDADO", xaxis = list(title = "DURACIÓN TRABAJO"), yaxis = list(title = "CANTIDAD"))
#Gráfica Trabajador extranjero
datos %>% count(default, foreign_worker) %>% plot_ly(x = ~default, y = ~n, color = ~foreign_worker, type = "bar") %>% layout(title = "TRABAJO EXTRANJERO VS CREDITO SALDADO", xaxis = list(title = "TRABAJO EXTRANJERO"), yaxis = list(title = "CANTIDAD"))
#Gráfica Proposito
datos %>% count(default, purpose) %>% plot_ly(x = ~default, y = ~n, color = ~purpose, type = "bar") %>% layout(title = "PROPOSITO CREDITO VS CREDITO SALDADO", xaxis = list(title = "PROPOSITO CREDITO"), yaxis = list(title = "CANTIDAD"))
Se puede observar que la variable edad, al tipo de trabajo de la persona y si trabaja en el extranjero no son influyentes en si se paga o no el crédito, ya que no hay mucha diferencia.
Respecto a la duración del trabajo de las personas, los datos suelen estar proporcionados, a excepción de los que tienen un trabajo de entre 1 - 4 años, que suelen ser menos morosos.
Finalmente, observamos que la gente suele pedir crédito para comprar un coche nuevo, muebles y radio/tv.
Ahora se va a comparar campos referentes al tipo de propiedades/inversiones
#Gráfica Antiguedad Residencia
plot_ly(datos, y = ~residence_history, color = ~default, type = "box") %>% layout(title = "ANTIGUEDAD RESIDENCIA VS CREDITO SALDADO")
#Gráfica Tipo de Vivienda
datos %>% count(default, housing) %>% plot_ly(x = ~default, y = ~n, color = ~housing, type = "bar") %>% layout(title = "TIPO DE VIVIENDA VS CREDITO SALDADO", xaxis = list(title = "TIPO DE VIVIENDA"), yaxis = list(title = "CANTIDAD"))
#Gráfica Propiedades o seguros
datos %>% count(default, property) %>% plot_ly(x = ~default, y = ~n, color = ~property, type = "bar") %>% layout(title = "PROPIEDADES VS CREDITO SALDADO", xaxis = list(title = "PROPIEDADES"), yaxis = list(title = "CANTIDAD"))
#Gráfica Otros planes de pago
datos %>% count(default, installment_plan) %>% plot_ly(x = ~default, y = ~n, color = ~installment_plan, type = "bar") %>% layout(title = "OTROS PLANES DE PAGO VS CREDITO SALDADO", xaxis = list(title = "OTROS PLANES DE PAGO"), yaxis = list(title = "CANTIDAD"))
Vemos que la antigüedad de la residencia y que tengan otras formas de pago no son unos factores que influyentes en la morosidad, moviéndose el grueso de estos valores entre 2,4 en la antigüedad de la residencia y que no dispongan de otro planes de pago.
Por otro lado, el tipo de propiedad que tienen son mayoritariamente en propiedad, no existiendo (a simple vista) ninguna variación rara respecto a la morosidad.
Finalmente, las personas que pagan el crédito suelen tener en mayor o menor medida alguna otra propiedad, mientras los que o tienen otro tipo de propiedad o del tipo “real estate” suelen ser los menos morosos.
Ahora se va a comparar como funcionan el resto de las variables.
#Gráfica Historial de credito
datos %>% count(default, credit_history) %>% plot_ly(x = ~default, y = ~n, color = ~credit_history, type = "bar") %>% layout(title = "HISTORIAL DE CREDITO VS CREDITO SALDADO", xaxis = list(title = "HISTORIAL DE CREDITO"), yaxis = list(title = "CANTIDAD"))
#Gráfica Otros deudores
datos %>% count(default, other_debtors) %>% plot_ly(x = ~default, y = ~n, color = ~other_debtors, type = "bar") %>% layout(title = "OTROS DEUDORES VS CREDITO SALDADO", xaxis = list(title = "OTROS DEUDORES"), yaxis = list(title = "CANTIDAD"))
#Gráfica Tasa de cuotas
plot_ly(datos, y = ~installment_rate, color = ~default, type = "box") %>% layout(title = "TASA DE CUOTAS VS CREDITO SALDADO")
#Gráfica Número de créditos
plot_ly(datos, y = ~existing_credits, color = ~default, type = "box") %>% layout(title = "NUMERO DE CREDITOS VS CREDITO SALDADO")
#Gráfica Telefono
datos %>% count(default, telephone) %>% plot_ly(x = ~default, y = ~n, color = ~telephone, type = "bar") %>% layout(title = "TELEFONO VS CREDITO SALDADO", xaxis = list(title = "TELEFONO"), yaxis = list(title = "CANTIDAD"))
Respecto a al campo historial de créditos, podemos ver que mayormente tienen algún otro crédito pagado anteriormente, sin embargo, los no morosos tienen algún otro crédito en estado crítico. En relación a este campo esta el numero de créditos, que aunque no es influyente, las personas suelen tener entre 1 y 2 créditos a la vez.
La tasa de cuotas de los créditos no algo influyente en la morosidad, pero suele ir mayormente de entre los 2 a 4 años, otro campo en las misma circunstancia es el poseer numero de teléfono, ya que hay mas casos en los que no se poseen.
Ahora, para comprobar los campos mas y menos influyente de una manera numérica, se hará unas pruebas estadísticas de significancia, para así determinar si se puede descartar algún campo. Para ellos se mirarán las proporciones, y luego se calculará los coeficientes V de Cramér y Phi.
if(!require(DescTools)){
install.packages('DescTools', repos='http://cran.us.r-project.org')
library(DescTools)
}
#Campo checking_balance
tabla_aux <- table(datos$checking_balance,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## < 0 DM 0.4927007 0.5072993
## > 200 DM 0.2222222 0.7777778
## 1 - 200 DM 0.3903346 0.6096654
## unknown 0.1167513 0.8832487
Phi(tabla_aux)
## [1] 0.3517399
CramerV(tabla_aux)
## [1] 0.3517399
El tipo de asociación es media, por lo que se deja el campo.
#Campo months_loan_duration
tabla_aux <- table(datos$months_loan_duration,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 4 0.0000000 1.0000000
## 5 0.0000000 1.0000000
## 6 0.1200000 0.8800000
## 7 0.0000000 1.0000000
## 8 0.1428571 0.8571429
## 9 0.2857143 0.7142857
## 10 0.1071429 0.8928571
## 11 0.0000000 1.0000000
## 12 0.2737430 0.7262570
## 13 0.0000000 1.0000000
## 14 0.2500000 0.7500000
## 15 0.1875000 0.8125000
## 16 0.5000000 0.5000000
## 18 0.3716814 0.6283186
## 20 0.1250000 0.8750000
## 21 0.3000000 0.7000000
## 22 0.0000000 1.0000000
## 24 0.3043478 0.6956522
## 26 0.0000000 1.0000000
## 27 0.3846154 0.6153846
## 28 0.3333333 0.6666667
## 30 0.3250000 0.6750000
## 33 0.3333333 0.6666667
## 36 0.4457831 0.5542169
## 39 0.2000000 0.8000000
## 40 1.0000000 0.0000000
## 42 0.2727273 0.7272727
## 45 0.8000000 0.2000000
## 47 0.0000000 1.0000000
## 48 0.5833333 0.4166667
## 54 0.5000000 0.5000000
## 60 0.4615385 0.5384615
## 72 1.0000000 0.0000000
Phi(tabla_aux)
## [1] 0.2808682
CramerV(tabla_aux)
## [1] 0.2808682
El tipo de asociación es media/baja, por lo que se deja el campo.
#Campo credit_history
tabla_aux <- table(datos$credit_history,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## critical 0.1706485 0.8293515
## delayed 0.3181818 0.6818182
## fully repaid 0.6250000 0.3750000
## fully repaid this bank 0.5714286 0.4285714
## repaid 0.3188679 0.6811321
Phi(tabla_aux)
## [1] 0.2483775
CramerV(tabla_aux)
## [1] 0.2483775
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo purpose
tabla_aux <- table(datos$purpose,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## business 0.3505155 0.6494845
## car (new) 0.3803419 0.6196581
## car (used) 0.1650485 0.8349515
## domestic appliances 0.3333333 0.6666667
## education 0.4400000 0.5600000
## furniture 0.3204420 0.6795580
## others 0.4166667 0.5833333
## radio/tv 0.2214286 0.7785714
## repairs 0.3636364 0.6363636
## retraining 0.1111111 0.8888889
Phi(tabla_aux)
## [1] 0.1826375
CramerV(tabla_aux)
## [1] 0.1826375
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo amount
tabla_aux <- table(datos$amount,datos$default )
Phi(tabla_aux)
## [1] 0.9652699
CramerV(tabla_aux)
## [1] 0.9652699
El tipo de asociación es muy alta, por lo que se deja el campo.
#Campo savings_balance
tabla_aux <- table(datos$savings_balance,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## < 100 DM 0.3598673 0.6401327
## > 1000 DM 0.1250000 0.8750000
## 101 - 500 DM 0.3300971 0.6699029
## 501 - 1000 DM 0.1746032 0.8253968
## unknown 0.1748634 0.8251366
Phi(tabla_aux)
## [1] 0.1899972
CramerV(tabla_aux)
## [1] 0.1899972
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo employment_length
tabla_aux <- table(datos$employment_length,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## > 7 yrs 0.2529644 0.7470356
## 0 - 1 yrs 0.4069767 0.5930233
## 1 - 4 yrs 0.3067847 0.6932153
## 4 - 7 yrs 0.2241379 0.7758621
## unemployed 0.3709677 0.6290323
Phi(tabla_aux)
## [1] 0.1355296
CramerV(tabla_aux)
## [1] 0.1355296
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo installment_rate
tabla_aux <- table(datos$installment_rate,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 1 0.2500000 0.7500000
## 2 0.2683983 0.7316017
## 3 0.2866242 0.7133758
## 4 0.3340336 0.6659664
Phi(tabla_aux)
## [1] 0.07400535
CramerV(tabla_aux)
## [1] 0.07400535
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo personal_status
tabla_aux <- table(datos$personal_status,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## divorced male 0.4000000 0.6000000
## female 0.3516129 0.6483871
## married male 0.2717391 0.7282609
## single male 0.2664234 0.7335766
Phi(tabla_aux)
## [1] 0.09800619
CramerV(tabla_aux)
## [1] 0.09800619
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo other_debtors
tabla_aux <- table(datos$other_debtors,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## co-applicant 0.4390244 0.5609756
## guarantor 0.1923077 0.8076923
## none 0.2998897 0.7001103
Phi(tabla_aux)
## [1] 0.08151912
CramerV(tabla_aux)
## [1] 0.08151912
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo residence_history
tabla_aux <- table(datos$residence_history,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 1 0.2769231 0.7230769
## 2 0.3149351 0.6850649
## 3 0.2885906 0.7114094
## 4 0.3002421 0.6997579
Phi(tabla_aux)
## [1] 0.02737328
CramerV(tabla_aux)
## [1] 0.02737328
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo property
tabla_aux <- table(datos$property,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## building society savings 0.3060345 0.6939655
## other 0.3072289 0.6927711
## real estate 0.2127660 0.7872340
## unknown/none 0.4350649 0.5649351
Phi(tabla_aux)
## [1] 0.1540115
CramerV(tabla_aux)
## [1] 0.1540115
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo age
tabla_aux <- table(datos$age,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 19 0.50000000 0.50000000
## 20 0.35714286 0.64285714
## 21 0.35714286 0.64285714
## 22 0.40740741 0.59259259
## 23 0.41666667 0.58333333
## 24 0.43181818 0.56818182
## 25 0.46341463 0.53658537
## 26 0.28000000 0.72000000
## 27 0.25490196 0.74509804
## 28 0.34883721 0.65116279
## 29 0.40540541 0.59459459
## 30 0.27500000 0.72500000
## 31 0.28947368 0.71052632
## 32 0.26470588 0.73529412
## 33 0.39393939 0.60606061
## 34 0.34375000 0.65625000
## 35 0.15000000 0.85000000
## 36 0.15384615 0.84615385
## 37 0.27586207 0.72413793
## 38 0.16666667 0.83333333
## 39 0.28571429 0.71428571
## 40 0.24000000 0.76000000
## 41 0.23529412 0.76470588
## 42 0.36363636 0.63636364
## 43 0.29411765 0.70588235
## 44 0.29411765 0.70588235
## 45 0.20000000 0.80000000
## 46 0.22222222 0.77777778
## 47 0.29411765 0.70588235
## 48 0.25000000 0.75000000
## 49 0.07142857 0.92857143
## 50 0.25000000 0.75000000
## 51 0.12500000 0.87500000
## 52 0.11111111 0.88888889
## 53 0.71428571 0.28571429
## 54 0.20000000 0.80000000
## 55 0.37500000 0.62500000
## 56 0.00000000 1.00000000
## 57 0.33333333 0.66666667
## 58 0.40000000 0.60000000
## 59 0.33333333 0.66666667
## 60 0.50000000 0.50000000
## 61 0.42857143 0.57142857
## 62 0.00000000 1.00000000
## 63 0.12500000 0.87500000
## 64 0.00000000 1.00000000
## 65 0.20000000 0.80000000
## 66 0.40000000 0.60000000
## 67 0.00000000 1.00000000
## 68 0.66666667 0.33333333
## 70 0.00000000 1.00000000
## 74 0.25000000 0.75000000
## 75 0.00000000 1.00000000
Phi(tabla_aux)
## [1] 0.2397444
CramerV(tabla_aux)
## [1] 0.2397444
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo installment_plan
tabla_aux <- table(datos$installment_plan,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## bank 0.4100719 0.5899281
## none 0.2751843 0.7248157
## stores 0.4042553 0.5957447
Phi(tabla_aux)
## [1] 0.1133101
CramerV(tabla_aux)
## [1] 0.1133101
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo housing
tabla_aux <- table(datos$housing,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## for free 0.4074074 0.5925926
## own 0.2608696 0.7391304
## rent 0.3910615 0.6089385
Phi(tabla_aux)
## [1] 0.1349068
CramerV(tabla_aux)
## [1] 0.1349068
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo existing_credits
tabla_aux <- table(datos$existing_credits,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 1 0.3159558 0.6840442
## 2 0.2762763 0.7237237
## 3 0.2142857 0.7857143
## 4 0.3333333 0.6666667
Phi(tabla_aux)
## [1] 0.05168364
CramerV(tabla_aux)
## [1] 0.05168364
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo dependents
tabla_aux <- table(datos$dependents,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## 1 0.3005917 0.6994083
## 2 0.2967742 0.7032258
Phi(tabla_aux)
## [1] 0.003014853
CramerV(tabla_aux)
## [1] 0.003014853
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo telephone
tabla_aux <- table(datos$telephone,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## none 0.3137584 0.6862416
## yes 0.2797030 0.7202970
Phi(tabla_aux)
## [1] 0.03646619
CramerV(tabla_aux)
## [1] 0.03646619
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo foreign_worker
tabla_aux <- table(datos$foreign_worker,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## no 0.1081081 0.8918919
## yes 0.3073728 0.6926272
Phi(tabla_aux)
## [1] 0.0820795
CramerV(tabla_aux)
## [1] 0.0820795
El tipo de asociación es baja, por lo que se descarta el campo.
#Campo job
tabla_aux <- table(datos$job,datos$default )
prop.table(tabla_aux, margin = 1)
##
## MOROSO NO MOROSO
## mangement self-employed 0.3445946 0.6554054
## skilled employee 0.2952381 0.7047619
## unemployed non-resident 0.3181818 0.6818182
## unskilled resident 0.2800000 0.7200000
Phi(tabla_aux)
## [1] 0.04341838
CramerV(tabla_aux)
## [1] 0.04341838
El tipo de asociación es baja, por lo que se descarta el campo.
Siguiendo los valores de V de Cramer y Phi, los valores entre 0.1 y 0.3 nos indican que la asociación estadística es baja, y entre 0.3 y 0.5 se puede considerar una asociación media. Finalmente, si los valores fueran superiores a 0.5, la asociación estadística entre las variables sería alta.
#Se hace una copia de los datos
datos_CRAMER_PHI_ALTO <- datos
#Se elimina los campos
datos_CRAMER_PHI_ALTO$dependents <- NULL
datos_CRAMER_PHI_ALTO$telephone <- NULL
datos_CRAMER_PHI_ALTO$job <- NULL
datos_CRAMER_PHI_ALTO$foreign_worker <- NULL
datos_CRAMER_PHI_ALTO$housing <- NULL
datos_CRAMER_PHI_ALTO$existing_credits <- NULL
datos_CRAMER_PHI_ALTO$installment_plan <- NULL
datos_CRAMER_PHI_ALTO$property <- NULL
datos_CRAMER_PHI_ALTO$purpose <- NULL
datos_CRAMER_PHI_ALTO$savings_balance <- NULL
datos_CRAMER_PHI_ALTO$employment_length <- NULL
datos_CRAMER_PHI_ALTO$installment_rate <- NULL
datos_CRAMER_PHI_ALTO$personal_status <- NULL
datos_CRAMER_PHI_ALTO$other_debtors <- NULL
datos_CRAMER_PHI_ALTO$residence_history <- NULL
datos_CRAMER_PHI_ALTO$age <- NULL
datos_CRAMER_PHI_ALTO$credit_history <- NULL
Ahora para proceder a preparar los datos, la primera cosa que debemos hacer es desordenar los datos.
set.seed(1)
data_random <- datos_CRAMER_PHI_ALTO[sample(nrow(datos_CRAMER_PHI_ALTO)),]
Como debemos dividir el conjunto de datos en dos grupos: entrenamiento y test, y al no existir un conjunto complementario ni proporción fijada, se hará 2/3 de los datos para el entrenamiento y 1/3 de los datos para el test.
La variable por la que clasificaremos es el campo de si la persona ha pagado el crédito o no, que está en la decimoséptima columna. De esta forma, tendremos un conjunto de datos para el entrenamiento y uno para la validación
set.seed(666)
y <- data_random[,4]
X <- data_random
X[,4] <- NULL
De forma dinámica podemos definir una forma de separar los datos en función de un parámetro, en este caso del “split_prop”. Definimos un parámetro que controla el split de forma dinámica en el test.
split_prop <- 3
max_split<-floor(nrow(X)/split_prop)
tr_limit <- nrow(X)-max_split
ts_limit <- nrow(X)-max_split+1
trainX <- X[1:tr_limit,]
trainy <- y[1:tr_limit]
testX <- X[(ts_limit+1):nrow(X),]
testy <- y[(ts_limit+1):nrow(X)]
En la segunda opción podemos crear directamente un rango utilizando el mismo parámetro anterior.
split_prop <- 3
indexes = sample(1:nrow(datos), size=floor(((split_prop-1)/split_prop)*nrow(datos)))
trainX<-X[indexes,]
trainy<-y[indexes]
testX<-X[-indexes,]
testy<-y[-indexes]
Al extraer aleatoriamente los datos, se hará un análisis mínimo de los datos para asegurarnos de no obtener clasificadores sesgados por los valores que contiene cada muestra.
En este caso, verificaremos que la proporción de morosos es más o menos constante en los dos conjuntos.
summary(trainX);
## checking_balance months_loan_duration amount
## Length:666 Min. : 4.00 Min. : 250
## Class :character 1st Qu.:12.00 1st Qu.: 1389
## Mode :character Median :18.00 Median : 2324
## Mean :21.05 Mean : 3280
## 3rd Qu.:24.00 3rd Qu.: 3973
## Max. :60.00 Max. :18424
summary(trainy)
## MOROSO NO MOROSO
## 204 462
summary(testX)
## checking_balance months_loan_duration amount
## Length:334 Min. : 4.00 Min. : 276
## Class :character 1st Qu.:12.00 1st Qu.: 1303
## Mode :character Median :18.00 Median : 2274
## Mean :20.61 Mean : 3254
## 3rd Qu.:24.00 3rd Qu.: 3970
## Max. :72.00 Max. :15945
summary(testy)
## MOROSO NO MOROSO
## 96 238
Se puede verificar, que hay aproximadamente la misma proporción en el conjunto de entrenamiento y de test. Siendo “NO” un poco mas del doble que “SI” en ambos casos.
Ya que tenemos los conjuntos preparados, se crea el árbol de decisión con los datos de entrenamiento.
trainy = as.factor(trainy)
model <- C50::C5.0(trainX, trainy,rules=TRUE )
summary(model)
##
## Call:
## C5.0.default(x = trainX, y = trainy, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Sat Dec 11 10:50:06 2021
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 666 cases (4 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (37/6, lift 2.7)
## checking_balance = < 0 DM
## months_loan_duration > 11
## amount <= 1391
## -> class MOROSO [0.821]
##
## Rule 2: (24/4, lift 2.6)
## checking_balance in {1 - 200 DM, < 0 DM}
## amount > 8588
## -> class MOROSO [0.808]
##
## Rule 3: (34/10, lift 2.3)
## checking_balance = < 0 DM
## months_loan_duration > 30
## -> class MOROSO [0.694]
##
## Rule 4: (160/72, lift 1.8)
## checking_balance = < 0 DM
## months_loan_duration > 11
## -> class MOROSO [0.549]
##
## Rule 5: (298/39, lift 1.2)
## checking_balance in {unknown, > 200 DM}
## -> class NO MOROSO [0.867]
##
## Rule 6: (629/182, lift 1.0)
## amount <= 8588
## -> class NO MOROSO [0.710]
##
## Default class: NO MOROSO
##
##
## Evaluation on training data (666 cases):
##
## Rules
## ----------------
## No Errors
##
## 6 152(22.8%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 71 133 (a): class MOROSO
## 19 443 (b): class NO MOROSO
##
##
## Attribute usage:
##
## 98.05% amount
## 71.62% checking_balance
## 24.02% months_loan_duration
##
##
## Time: 0.0 secs
El árbol obtenido clasifica erróneamente 152 de los 666 casos dados, una tasa de error del 22.8%.
A partir del árbol de decisión, se pueden extraer las siguientes 6 reglas de:
checking_balance = < 0 DM && months_loan_duration > 11 && amount <= 1391 –> MOROSO. Validez: 82,1%
checking_balance entre los valores {1 - 200 DM, < 0 DM} && amount > 8588 –> MOROSO. Validez: 80,8%
checking_balance = < 0 DM && months_loan_duration > 30 –> MOROSO. Validez: 69,4%
checking_balance = < 0 DM && months_loan_duration > 11 –> MOROSO. Validez: 54,9%
checking_balance entre {unknown, > 200 DM} –> NO MOROSO. Validez: 86,7%
amount <= 8588 –> NO MOROSO. Validez: 71%
A continuación, mostramos el árbol obtenido.
model <- C50::C5.0(trainX, trainy)
plot(model)
Una vez tenemos el modelo, podemos comprobar su calidad prediciendo la clase para los datos de prueba que nos hemos reservado al principio.
predicted_model <- predict( model, testX, type="class" )
print(sprintf("La precisión del árbol es: %.4f %%",100*sum(predicted_model == testy) / length(predicted_model)))
## [1] "La precisión del árbol es: 73.0539 %"
Cuando hay pocas clases, la calidad de la predicción se puede analizar mediante una matriz de confusión que identifica los tipos de errores cometidos.
mat_conf<-table(testy,Predicted=predicted_model)
mat_conf
## Predicted
## testy MOROSO NO MOROSO
## MOROSO 25 71
## NO MOROSO 19 219
Para tener información más completa se usará el paquete gmodels.
if(!require(gmodels)){
install.packages('gmodels', repos='http://cran.us.r-project.org')
library(gmodels)
}
CrossTable(testy, predicted_model,prop.chisq = FALSE, prop.c = FALSE, prop.r =FALSE,dnn = c('Reality', 'Prediction'))
##
##
## Cell Contents
## |-------------------------|
## | N |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 334
##
##
## | Prediction
## Reality | MOROSO | NO MOROSO | Row Total |
## -------------|-----------|-----------|-----------|
## MOROSO | 25 | 71 | 96 |
## | 0.075 | 0.213 | |
## -------------|-----------|-----------|-----------|
## NO MOROSO | 19 | 219 | 238 |
## | 0.057 | 0.656 | |
## -------------|-----------|-----------|-----------|
## Column Total | 44 | 290 | 334 |
## -------------|-----------|-----------|-----------|
##
##
Como se puede observar, el árbol tiene una precisión de un poco mas del 73%, lo que está bastante bien. Lo que voy a comprobar ahora es la precisión del árbol con todas las variables, ya que se ha descartado la inmensa mayoría.
#Asignamos los datos
set.seed(1)
data_random_completos <- datos[sample(nrow(datos)),]
#Separamos los valores
set.seed(666)
y_completo <- data_random_completos[,17]
X_completo <- data_random_completos
X_completo[,17] <- NULL
#Separamos los campos
split_prop <- 3
max_split<-floor(nrow(X_completo)/split_prop)
tr_limit <- nrow(X_completo)-max_split
ts_limit <- nrow(X_completo)-max_split+1
trainX <- X_completo[1:tr_limit,]
trainy <- y_completo[1:tr_limit]
testX <- X_completo[(ts_limit+1):nrow(X_completo),]
testy <- y_completo[(ts_limit+1):nrow(X_completo)]
split_prop <- 3
indexes = sample(1:nrow(datos), size=floor(((split_prop-1)/split_prop)*nrow(datos)))
trainX<-X_completo[indexes,]
trainy<-y_completo[indexes]
testX<-X_completo[-indexes,]
testy<-y_completo[-indexes]
#Se crea el arbol de decisión
trainy = as.factor(trainy)
model <- C50::C5.0(trainX, trainy,rules=TRUE )
#Se obtiene la precision del arbol
predicted_model <- predict( model, testX, type="class" )
print(sprintf("La precisión del árbol con todos los campos es: %.4f %%",100*sum(predicted_model == testy) / length(predicted_model)))
## [1] "La precisión del árbol con todos los campos es: 73.6527 %"
La variación de precisión es mínima con todos los campos respecto al modelo con campos simplificados.
Se puede concluir que la capacidad de predicción del árbol es bastante buena, y que como se ha comprobado un análisis inicial de los campos, pueden ayudar a simplificar mucho la creación del árbol.
Ahora se va a implementar un modelo complementario con el paquete rpart de R. Lo primero es crear nuestro conjunto de entrenamiento y de prueba (Ejemplo de arbol basado en: https://rpubs.com/jboscomendoza/arboles_decision_clasificacion)
library(tidyverse)
library(rpart)
library(rpart.plot)
library(caret)
#Separamos los campos
set.seed(666)
y <- data_random[,4]
X <- data_random
X[,4] <- NULL
#Separamos los valores
split_prop <- 3
max_split<-floor(nrow(X)/split_prop)
tr_limit <- nrow(X)-max_split
ts_limit <- nrow(X)-max_split+1
trainX <- X[1:tr_limit,]
trainy <- y[1:tr_limit]
testX <- X[(ts_limit+1):nrow(X),]
testy <- y[(ts_limit+1):nrow(X)]
split_prop <- 3
indexes = sample(1:nrow(datos), size=floor(((split_prop-1)/split_prop)*nrow(datos)))
trainX<-X[indexes,]
trainy<-y[indexes]
testX<-X[-indexes,]
testy<-y[-indexes]
Usamos la función rpart de rpart para entrenar nuestro modelo. Esta función nos pide una fórmula para especificar la variable objetivo de la clasificación. La fórmula que usaremos es tipo ~ ., la cual expresa que intentaremos clasificar tipo usando a todas las demás variables como predictoras.
arbol <- rpart(formula = trainy ~ ., data = trainX)
arbol
## n= 666
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 666 204 NO MOROSO (0.30630631 0.69369369)
## 2) checking_balance=< 0 DM,1 - 200 DM 368 165 NO MOROSO (0.44836957 0.55163043)
## 4) months_loan_duration>=22.5 167 72 MOROSO (0.56886228 0.43113772)
## 8) amount< 1549.5 13 1 MOROSO (0.92307692 0.07692308) *
## 9) amount>=1549.5 154 71 MOROSO (0.53896104 0.46103896)
## 18) amount>=8760.5 21 4 MOROSO (0.80952381 0.19047619) *
## 19) amount< 8760.5 133 66 NO MOROSO (0.49624060 0.50375940)
## 38) months_loan_duration>=26.5 83 36 MOROSO (0.56626506 0.43373494)
## 76) amount< 2165.5 7 1 MOROSO (0.85714286 0.14285714) *
## 77) amount>=2165.5 76 35 MOROSO (0.53947368 0.46052632)
## 154) amount>=5981.5 25 8 MOROSO (0.68000000 0.32000000) *
## 155) amount< 5981.5 51 24 NO MOROSO (0.47058824 0.52941176)
## 310) months_loan_duration>=37.5 12 4 MOROSO (0.66666667 0.33333333) *
## 311) months_loan_duration< 37.5 39 16 NO MOROSO (0.41025641 0.58974359) *
## 39) months_loan_duration< 26.5 50 19 NO MOROSO (0.38000000 0.62000000) *
## 5) months_loan_duration< 22.5 201 70 NO MOROSO (0.34825871 0.65174129)
## 10) amount< 1285 71 34 NO MOROSO (0.47887324 0.52112676)
## 20) months_loan_duration>=11 47 17 MOROSO (0.63829787 0.36170213) *
## 21) months_loan_duration< 11 24 4 NO MOROSO (0.16666667 0.83333333) *
## 11) amount>=1285 130 36 NO MOROSO (0.27692308 0.72307692)
## 22) amount>=6138.5 13 5 MOROSO (0.61538462 0.38461538) *
## 23) amount< 6138.5 117 28 NO MOROSO (0.23931624 0.76068376) *
## 3) checking_balance=> 200 DM,unknown 298 39 NO MOROSO (0.13087248 0.86912752) *
Lo anterior muestra el esquema de nuestro árbol de clasificación. Cada inciso nos indica un nodo y la regla de clasificación que le corresponde. Siguiendo estos nodos, podemos llegar a las hojas del árbol, que corresponde a la clasificación de nuestros datos.
Todo lo anterior resulta mucho más claro si lo visualizamos.
rpart.plot(arbol)
Ahora generar un vector con los valores predichos por el modelo que hemos entrenado y mostramos la matriz de confusión.
prediccion <- predict(arbol, newdata = testX, type = "class")
confusionMatrix(prediccion, testy)
## Confusion Matrix and Statistics
##
## Reference
## Prediction MOROSO NO MOROSO
## MOROSO 34 33
## NO MOROSO 62 205
##
## Accuracy : 0.7156
## 95% CI : (0.6639, 0.7633)
## No Information Rate : 0.7126
## P-Value [Acc > NIR] : 0.479311
##
## Kappa : 0.2369
##
## Mcnemar's Test P-Value : 0.004069
##
## Sensitivity : 0.3542
## Specificity : 0.8613
## Pos Pred Value : 0.5075
## Neg Pred Value : 0.7678
## Prevalence : 0.2874
## Detection Rate : 0.1018
## Detection Prevalence : 0.2006
## Balanced Accuracy : 0.6078
##
## 'Positive' Class : MOROSO
##
Y calculamos la precisión del árbol generado:
predicted_model <- predict( arbol, testX, type="class" )
print(sprintf("La precisión del árbol es: %.4f %%",100*sum(predicted_model == testy) / length(predicted_model)))
## [1] "La precisión del árbol es: 71.5569 %"
Se observa de manera clara, que las precisiones de los arboles generados no suele variar mucho. No obstante, se va a hacer otro ultimo árbol basado en el ejemplo del punto 2.5 del enunciado que incorpora como novedad “adaptative boosting”.
modelo2 <- C50::C5.0(trainX, trainy, trials = 10)
plot(modelo2)
Vemos a continuación cómo son las predicciones del nuevo árbol:
predicted_model2 <- predict( modelo2, testX, type="class" )
print(sprintf("La precisión del árbol es: %.4f %%",100*sum(predicted_model2 == testy) / length(predicted_model2)))
## [1] "La precisión del árbol es: 72.4551 %"
Observamos como se modifica levemente la precisión del modelo a peor.
Como conclusión de esta PEC3, me he dado cuenta la importancia de realizar un análisis previo de las variables y su relación con el propósito que queremos buscar, ya que no es lo mismo trabajar con 21 campos que con 5 (como ha ocurrido en este caso), y aunque el modelo de predicción con menos campos ha sido muy ligeramente inferior (menor que un 1%) la forma de comprender las reglas y visualizar el árbol ha sido mucho más fácil.
Respecto a los tres métodos de arboles (dos realmente ya que uno era con una variación), creo que pueden ser complementarios, ya que uno aporta un enfoque mas simples respecto a las reglas, y otro una visualización más clara.